Apply Models

Generate spatial predictions across a landscape

Hong Jhun Sim

2022-09-05

This article outlines the steps to predict the number of species (species richness) for a particular animal group, across a target landscape. The pixel-based predictions are then visualised in the form of a spatial heatmap. Begin by loading the required packages to run the analysis:

library("biodivercity")
library("dplyr") # to process/wrangle data
library("tmap") # for visualisation

The biodivercity package contains predictive models (lme4::glmer()) and data pre-processing workflow ‘recipes’ (recipes::recipe()) that were built for each of four animal groups surveyed in Singapore. Refer to vignette("build-models") for more details on these models, which were built using the example data in this package. In this article, the model/recipe used to predict the number of bird species will be used as an example. Alternatively, the user may provide their own model/recipe for the analysis.

Load the model/recipe objects built exclusively from manually-generated landscape data:

filepath <- system.file("extdata", "models-manually-mapped.Rdata", package = "biodivercity")
load(filepath)

Landscape data within the target area of interest will be used to make spatial predictions. They will be summarised according to the predictor variables present in the models. For example, predictor variables in the bird model can be extracted as follows. Note that the prefix r<value>m denotes the buffer radius that the particular landscape data was summarised within, and _man_ denotes that the variables were generated from manually-generated data.

predictors_birds <- models_birds %>% # a list object
 lapply(function(x) names(x@frame)) %>% # extract variable names in model
  unlist() %>%
  unique() %>% 
  stringr::str_subset("(?<=^r)\\d+.*") # predictor variables start with the "r<value>m_"

predictors_birds

## [1] "r126m_man_natveg_pland"    "r50m_man_buildingFA_ratio"
## [3] "r50m_man_laneDensity"      "r50m_man_tree_sprich"     
## [5] "r50m_man_water_pland"      "r50m_man_shrub_pland"     
## [7] "r50m_man_turf_pland"       "r50m_man_shrub_sprich"

The following list describes the manually-generated landscape components used to build models provided in this package, their respective vector format, as well as all the possible predictor variables that may be summarised:

• Natural vegetation (polygons)
  • Percentage of landscape area (natveg_pland)
• Trees (points) with species name (per point)
  • Species richness (tree_sprich)
• Shrubs (polygons) with species name (per polygon)
  • Percentage of landscape area (shrub_pland)
  • Species richness (shrub_sprich)
• Turf (polygons)
  • Percentage of landscape area (turf_pland)
• Water (polygons)
  • Percentage of landscape area (water_pland)
• Buildings (polygons) each with the number of levels
  • Floor area ratio (buildingFA_ratio)
    
  • Average number of levels (buildingAvgLvl)
• Roads (lines) each with the number of lanes
  • Lane density (laneDensity)

To demonstrate how spatial predictions may be generated, load the example landscape data from within the Punggol (PG) area in Singapore (Chong et al., 2014, 2019), visualised in the interactive map below. Note that there are no water bodies mapped in the target area.

filepath <- system.file("extdata", "pg_layers.Rdata", package = "biodivercity")
load(filepath)

Spatial predictions

To make spatial predictions, the target area is broken up into many smaller points (pixels), where landscape data will be summarised and predictions will be made. First, the function generate_grid() will be used to generate a grid over the target area (see Figure 1). The pixel resolution of the grid can be specified with the argument pixelsize_m.

Figure 1: Visualisation of the function generate_grid().

An optional argument innerbuffer_m is also provided. This limits output predictions to a smaller area within the target landscape (Figure 1). This is useful, for example, when the model requires broad-scale landscape data to make predictions, but landscape data that is available do not extend beyond the target area. Hence, to avoid inaccurate predictions that result from areas without landscape data, the distance value for this argument should correspond to the largest buffer radius present in the model variables. In this example, the largest radius for variables within the bird models is 126 metres:

# get the max radius among all predictor variables
max_radius <- predictors_birds %>% 
  stringr::str_extract("(?<=^r)\\d+") %>%  # extract values for buffer radii
  as.numeric() %>% 
  max(na.rm = TRUE)
  
max_radius

## [1] 126

With the max_radius defined, run the function generate_grid(). The output is a dataframe of points (see Figure 1), where landscape variables will be summarised within their respective distance buffers. In this example, we will generate a grid with a pixel resolution of 50 metres.

grid_points <- generate_grid(target_areas = bound, # target area
                             innerbuffer_m = max_radius, # exclude areas < 126 m from boundaries
                             pixelsize_m = 50) %>%  # pixel resolution
  rownames_to_column("point_id") # add unique identifier

grid_points # geometry column has been added

## Simple feature collection with 306 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 378068.1 ymin: 153518.8 xmax: 379468.1 ymax: 154968.8
## Projected CRS: WGS 84 / UTM zone 48N
## First 10 features:
##    point_id area                  geometry
## 1         1   PG POINT (379168.1 153518.8)
## 2         2   PG POINT (379118.1 153568.8)
## 3         3   PG POINT (379168.1 153568.8)
## 4         4   PG POINT (379068.1 153618.8)
## 5         5   PG POINT (379118.1 153618.8)
## 6         6   PG POINT (379168.1 153618.8)
## 7         7   PG POINT (379018.1 153668.8)
## 8         8   PG POINT (379068.1 153668.8)
## 9         9   PG POINT (379118.1 153668.8)
## 10       10   PG POINT (379168.1 153668.8)

Next, use the function calc_manual() to summarise each landscape component around each of the point locations in grid_points, and at specified buffer radii. For example, to use the bird models, we summarise the landscape data as follows. The predictor variables can be appended to the grid_points dataframe as new columns, corresponding to those present in predictors_birds.

predictors_buildings <- 
  calc_manual(vector = buildings, name = "buildings",
              points = grid_points, buffer_sizes = 50,
              building_levels = "levels") %>% 
   lapply(st_drop_geometry) %>% 
  bind_rows(.id = "radius_m") %>% 
  pivot_wider(names_from = "radius_m", 
              values_from = starts_with("man"),
              names_glue = "r{radius_m}m_{.value}")

predictors_roads <- 
  calc_manual(vector = roads, name = "roads",
              points = grid_points, buffer_sizes = 50,
              road_lanes = "lanes") %>% 
  lapply(st_drop_geometry) %>% 
  bind_rows(.id = "radius_m") %>% 
  pivot_wider(names_from = "radius_m", 
              values_from = starts_with("man"),
              names_glue = "r{radius_m}m_{.value}")

predictors_trees <- 
  calc_manual(vector = trees, name = "trees",
              points = grid_points, buffer_sizes = 50,
              plant_species = "species") %>% 
  lapply(st_drop_geometry) %>% 
  bind_rows(.id = "radius_m") %>% 
  pivot_wider(names_from = "radius_m", 
              values_from = starts_with("man"),
              names_glue = "r{radius_m}m_{.value}")

predictors_shrubs <- 
  calc_manual(vector = shrubs, name = "shrubs",
              points = grid_points, buffer_sizes = 50,
              plant_species = "species") %>% 
  lapply(st_drop_geometry) %>% 
  bind_rows(.id = "radius_m") %>% 
  pivot_wider(names_from = "radius_m", 
              values_from = starts_with("man"),
              names_glue = "r{radius_m}m_{.value}")

predictors_turf <- 
  calc_manual(vector = turf, name = "turf", 
              points = grid_points, buffer_sizes = 50) %>% 
  lapply(st_drop_geometry) %>% 
  bind_rows(.id = "radius_m") %>% 
  pivot_wider(names_from = "radius_m", 
              values_from = starts_with("man"),
              names_glue = "r{radius_m}m_{.value}")

predictors_natveg <- 
  calc_manual(vector = natveg, name = "natveg",
              points = grid_points, buffer_sizes = c(50, 126)) %>% 
  lapply(st_drop_geometry) %>% 
  bind_rows(.id = "radius_m") %>% 
  pivot_wider(names_from = "radius_m", 
              values_from = starts_with("man"),
              names_glue = "r{radius_m}m_{.value}")

predictors_water <- 
  calc_manual(vector = water, name = "water",
              points = grid_points, buffer_sizes = 50) %>% 
  lapply(st_drop_geometry) %>% 
  bind_rows(.id = "radius_m") %>% 
  pivot_wider(names_from = "radius_m", 
              values_from = starts_with("man"),
              names_glue = "r{radius_m}m_{.value}")


# combine all landscape predictors
grid_landscape <- grid_points %>% 
  inner_join(predictors_buildings) %>% 
  inner_join(predictors_roads) %>% 
  inner_join(predictors_trees) %>% 
  inner_join(predictors_shrubs) %>% 
  inner_join(predictors_turf) %>% 
  inner_join(predictors_natveg) %>% 
  inner_join(predictors_water)

grid_landscape

## Simple feature collection with 306 features and 16 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 378068.1 ymin: 153518.8 xmax: 379468.1 ymax: 154968.8
## Projected CRS: WGS 84 / UTM zone 48N
## First 10 features:
##    point_id area r50m_man_buildingArea_m2 r50m_man_buildingGFA_m2
## 1         1   PG                        0                       0
## 2         2   PG                        0                       0
## 3         3   PG                        0                       0
## 4         4   PG                        0                       0
## 5         5   PG                        0                       0
## 6         6   PG                        0                       0
## 7         7   PG                        0                       0
## 8         8   PG                        0                       0
## 9         9   PG                        0                       0
## 10       10   PG                        0                       0
##    r50m_man_buildingAvgLvl r50m_man_buildingFA_ratio r50m_man_laneLength_m
## 1                        0                         0                0.0000
## 2                        0                         0                0.0000
## 3                        0                         0                0.0000
## 4                        0                         0              207.2115
## 5                        0                         0                0.0000
## 6                        0                         0                0.0000
## 7                        0                         0              498.5300
## 8                        0                         0              582.2179
## 9                        0                         0              196.7592
## 10                       0                         0                0.0000
##    r50m_man_laneDensity r50m_man_tree_sprich r50m_man_tree_count
## 1            0.00000000                    0                   0
## 2            0.00000000                    0                   0
## 3            0.00000000                    0                   0
## 4            0.02638299                    1                   4
## 5            0.00000000                    0                   0
## 6            0.00000000                    0                   0
## 7            0.06347482                    1                   7
## 8            0.07413028                    2                  19
## 9            0.02505215                    2                   9
## 10           0.00000000                    0                   0
##    r50m_man_shrub_sprich r50m_man_shrub_pland r50m_man_turf_pland
## 1                      0       0.000000 [m^2]      0.000000 [m^2]
## 2                      0       0.000000 [m^2]      0.000000 [m^2]
## 3                      0       0.000000 [m^2]      0.000000 [m^2]
## 4                      1       2.442749 [m^2]      9.646878 [m^2]
## 5                      0       0.000000 [m^2]      0.000000 [m^2]
## 6                      0       0.000000 [m^2]      0.000000 [m^2]
## 7                      1       3.914387 [m^2]     62.940872 [m^2]
## 8                      1       6.544102 [m^2]     35.268007 [m^2]
## 9                      1       1.911862 [m^2]     13.373559 [m^2]
## 10                     0       0.000000 [m^2]      0.000000 [m^2]
##    r50m_man_natveg_pland r126m_man_natveg_pland r50m_man_water_pland
## 1       99.9543137 [m^2]         99.95431 [m^2]                    0
## 2       99.9543137 [m^2]         92.18319 [m^2]                    0
## 3       99.9543137 [m^2]         99.51171 [m^2]                    0
## 4       74.6022679 [m^2]         62.14158 [m^2]                    0
## 5       99.9519902 [m^2]         78.20742 [m^2]                    0
## 6       99.9543137 [m^2]         91.52636 [m^2]                    0
## 7        0.6084413 [m^2]         26.96329 [m^2]                    0
## 8       26.3762809 [m^2]         44.21388 [m^2]                    0
## 9       73.7248376 [m^2]         62.42091 [m^2]                    0
## 10      99.9543137 [m^2]         78.08538 [m^2]                    0
##                     geometry
## 1  POINT (379168.1 153518.8)
## 2  POINT (379118.1 153568.8)
## 3  POINT (379168.1 153568.8)
## 4  POINT (379068.1 153618.8)
## 5  POINT (379118.1 153618.8)
## 6  POINT (379168.1 153618.8)
## 7  POINT (379018.1 153668.8)
## 8  POINT (379068.1 153668.8)
## 9  POINT (379118.1 153668.8)
## 10 POINT (379168.1 153668.8)

Finally, use the function predict_heatmap() to predict the number of bird species (species richness) across the generated grid of spatial points, based on the predictor variables summarised in grid_landscape. The model object models_birds (suite of ‘best’ models) and recipe_birds (data pre-processing workflow) will be used to make the predictions at each point. Ensure that the argument pixelsize_m is similar to the value used in generate_grid().

bird_heatmap <- predict_heatmap(models = models_birds, 
                                recipe_data = recipe_birds,
                                points_topredict = grid_landscape, 
                                pixelsize_m = 50)

Visualisation

The continuous raster may be visualised as a heatmap:

tmap_mode("view")
tmap_options(max.raster = c(view = 1e8)) # increase max resolution to be visualised

tm_basemap(c("CartoDB.Positron", "OpenStreetMap")) +
  tm_shape(bird_heatmap, raster.downsample = FALSE) +
  tm_raster(title = "Number of bird species",
            style = "pretty",
            palette = "YlOrRd",
            alpha = 0.6) 

---

References

Chong KY, Teo S, Kurukulasuriya B, Chung YF, Giam X & Tan HTW (2019) The effects of landscape scale on greenery and traffic relationships with urban birds and butterflies. Urban Ecosystems, 22(5): 917–926.

Chong KY, Teo S, Kurukulasuriya B, Chung YF, Rajathurai S & Tan HTW (2014) Not all green is as good: Different effects of the natural and cultivated components of urban vegetation on bird and butterfly diversity. Biological Conservation, 171: 299–309.